---
title: "Negativity Biases and Political Ideology: A Comparative Test Across 17 Countries, Replication File"
output: html_notebook
---

The file that follows replicates all analyses in Patrick Fournier, Stuart Soroka and Lilach Nir (2020), Negativity Biases and Political Ideology: A Comparative Test Across 17 Countries, forthcoming in the American Political Science Review.  Analyses were done primiarly in STATA, with selected analyses for figures - and the figures themselves - drawn in R.  For this reason, we are distributing this RMarkdown file that pulls together all the various datasets and analyses.

The only graphic that is not produced using the distributed script is Figure 5, which offers descriptives for the IAPS pictures used in the study.  Those descriptives must be downloaded here: https://csea.phhp.ufl.edu.

There are five different datasets distributed with this RMarkdown file:

Sample-description.dta <- The respondent-level dataset including all responses, not just those with working physiological measures.

Individual-data.dta <- The working respondent-level dataset.

Stimulus-data.dta <- The respondent-stimulus-level panel dataset.

Time-series-data-photos.dta <- The times-series by-respondent panel dataset.

Time-series-data-videos.dta <- The times-series by-respondent panel dataset.

Analysis run in STATA are stored in STATA do-files with names that correspond to each of these datasets.  Those STATA do-files can be run in STATA, independent of this script.  Alternatively, the current file produces the figures alongside reading in and processing all of the STATA do-files using the RStata package.

## Setup

```{r setup}

rm(list=ls())

library(rio) ; library(plm) ; library(RStata) ; library(effects) ; library(knitr)
options("RStata.StataPath" = '/Applications/Stata/StataMP.app/Contents/MacOS/stata-mp')
options("RStata.StataVersion" = 15)

#setting countries
countries <- c("BR","CA.e","CA.f","CH","CN","DK","FR","GH","IN","IS.j","IS.p","IT","JP","NZ","RU","SE","SW","UK","US")
ncountries <- length(countries)

```

```{r datasets}

S1 <- import("Sample-description.dta") #respondent-level dataset
S1 <- S1[S1$country2!="",]
S <- import("Individual-data.dta") #respondent-level dataset
S <- S[S$country2!="",]
R <- import("Stimulus-data.dta") #respondent-level dataset
V <- import("Time-series-data-videos.dta") #respondent-level dataset
P <- import("Time-series-data-photos.dta") #respondent-level dataset

```

## Analysis of Sample-description.dta

```{r figure 1a}

TAB <- as.data.frame(table(S1$respnum))
TAB<- table(S1$female,S1$country2)
TAB <- TAB[1:2,ncountries:1]
{
pdf("figure.1a.pdf", width = 4, height = 6, bg="white") 
barplot(TAB, las=1, legend.text = c("Female", "Male"),args.legend = list(x = "right", bty = "n"), horiz=TRUE)
title(main="Respondent Sex")
title(xlab="Sample Size by Sex")
invisible(dev.off())
}

```

```{r figure 1b}

{
pdf("figure.1b.pdf", width = 4, height = 6, bg="white") 
plot(c(18,70),c(0,ncountries),ylim=c(.5,ncountries), ann=F,axes=F,type="n")
title(main="Respondent Age")
for (i in 1:ncountries) {
  Q <- quantile(S1$age[S1$country2==countries[i]],c(.25,.75),na.rm=T)
  arrows(Q[1],ncountries+1-i,Q[2],ncountries+1-i,code=3,length=.05,angle=90,col="gray")
}
for (i in 1:ncountries) {
  M <- mean(S1$age[S1$country2==countries[i]],na.rm=T) ; points(M,ncountries+1-i,pch=15,cex=1.2)
}
axis(1)
axis(2,at=c(1:ncountries),labels=rev(countries),las=1)
title(xlab="Age, in Years")
invisible(dev.off())
}

```

```{r figure 1c}

{
pdf("figure.1c.pdf", width = 4, height = 6, bg="white") 
plot(c(1,2,3,4,5,6,7),c(0,2,5,6,7,11,12),ylim=c(.5,ncountries), ann=F,axes=F,type="n")
title(main="English Proficiency")
for (i in 1:ncountries) {
  Q <- quantile(S1$english_prof[S1$country2==countries[i]],c(.25,.75),na.rm=T)
  arrows(Q[1],ncountries+1-i,Q[2],ncountries+1-i,code=3,length=.05,angle=90,col="gray")
}
for (i in 1:ncountries) {
  M <- mean(S1$english_prof[S1$country2==countries[i]]) ; points(M,ncountries+1-i,pch=15,cex=1.2)
}
axis(1)
axis(2,at=c(1:ncountries),labels=rev(countries),las=1)
title(xlab="English Proficiency")
invisible(dev.off())
}

```

```{r STATA sample-description analyses, as.is=T}

stata_src <- readLines("Sample-description.do")
stata(stata_src,data.in=S1)

```

## Analysis of Individual-data.dta

```{r Figure 2}

M <- S[,c("v_gslmc_11","v_gslmc_12","v_gslmc_13","v_gslmc_14")]
M <- rowMeans(M,na.rm=T)
S$mean.neg <- M

M <- S[,c("v_gslmc_15","v_gslmc_16","v_gslmc_17","v_gslmc_18")]
M <- rowMeans(M,na.rm=T)
S$mean.pos <- M

a1 <- S[,c("mean.neg","wp_right2_dichox")]
a1$tone <- "neg"
colnames(a1) <- c("gsl","lr","tone")
a2 <- S[,c("mean.pos","wp_right2_dichox")]
a2$tone <- "pos"
colnames(a2) <- c("gsl","lr","tone")
A <- rbind(a1,a2)
rm(a1,a2)
A$tone <- as.factor(A$tone)
A$lr <- as.factor(A$lr)

model1 <- lm(gsl ~ lr * tone, data=A)
eff1 <- effect("lr * tone",model1, typical=mean)
s <- cbind(eff1$x,fit=eff1$fit,lower=eff1$lower,higher=eff1$upper)
s <- s[order(s$lr,s$tone),]
s$range <- c(1,2,4,5)

{
pdf("figure2.pdf", width = 7, height = 5, bg="white") 
par(mar=c(5.1,6.1,1.1,2.1))
plot(s$range,s$fit,type="n",ann=F,axes=F,xlim=c(.5,5.5),ylim=c(-.16,.08))
arrows(s$range,s$low,s$range,s$high,angle=90,length=.05,code=3,col="gray",lwd=2,lty=1)
points(s$range,s$fit,pch=15,cex=1.5)
axis(1,col="white",at=c(1,2,4,5),labels=c("Negative\n Video","Positive\n Video","Negative\n Video", "Positive\n Video"))
axis(2,las=1)
mtext("Mean Normalized GSL",side=2,line=4)
mtext("Left-Leaning Participants              Right-Leaning Participants",side=1,line=3)
invisible(dev.off())
}

```

```{r Figure 6}

a1 <- S[,c("mean_ph_neg_all","wp_right2_dichox")]
a1$tone <- "neg"
colnames(a1) <- c("gsl","lr","tone")
a2 <- S[,c("mean_ph_pos_all","wp_right2_dichox")]
a2$tone <- "pos"
colnames(a2) <- c("gsl","lr","tone")
A <- rbind(a1,a2)
rm(a1,a2)
A$tone <- as.factor(A$tone)
A$lr <- as.factor(A$lr)

model1 <- lm(gsl ~ lr * tone, data=A)
eff1 <- effect("lr * tone",model1, typical=mean)
s <- cbind(eff1$x,fit=eff1$fit,lower=eff1$lower,higher=eff1$upper)
s <- s[order(s$lr,s$tone),]
s$range <- c(1,2,4,5)

{
pdf("figure6.pdf", width = 7, height = 5, bg="white") 
par(mar=c(5.1,6.1,1.1,2.1))
plot(s$range,s$fit,type="n",ann=F,axes=F,xlim=c(.5,5.5),ylim=c(-.05,.05))
arrows(s$range,s$low,s$range,s$high,angle=90,length=.05,code=3,col="gray",lwd=2,lty=1)
points(s$range,s$fit,pch=15,cex=1.5)
axis(1,col="white",at=c(1,2,4,5),labels=c("Negative\n Photos","Positive\n Photos","Negative\n Photos", "Positive\n Photos"))
axis(2,las=1)
mtext("Mean Normalized GSL",side=2,line=4)
mtext("Left-Leaning Participants              Right-Leaning Participants",side=1,line=3)
invisible(dev.off())
}

```

```{r Appendix Figure 1}

S$right <- NA
S$right[S$left_right2<.5] <- 0
S$right[S$left_right2==.5] <- .5
S$right[S$left_right2>.5] <- 1
tab <- RcmdrMisc::rowPercents(table(S$country2,S$right))
tab <- t(tab[,c(1:3)])
tab <- tab[,c(19:1)]

{
pdf("figureA1.pdf", width = 4, height = 6, bg="white") 
barplot(tab, las=1, legend.text = c("Left", "Center", "Right"),args.legend = list(x = "bottom", bty = "n", inset=c(0, -.23),ncol=3), horiz=TRUE)
mtext("% Participants by Ideological Category", side=1, line=2)
invisible(dev.off())
}

```

```{r STATA individual-data analyses, as.is=T}

stata_src <- readLines("Individual-data.do")
stata(stata_src,data.in=S)

```

## Analysis of Stimulus-data.dta

```{r Figure 3}

coefs <- as.data.frame(countries)
coefs$b <- NA
coefs$se <- NA
E <- R[R$stimtype=="single video" & R$local==0,]
E <- E[!is.na(E$vid_order),]
E$negativity <- E$v_negativity
Pd <- pdata.frame(E,index = c("resp", "vid_order"))
Pd$orderN <- as.numeric(Pd$vid_order)
for (i in countries){
  model <- plm(gslmc ~ negativity * wp_right2 + orderN + female + age + income + university,
               data=Pd[Pd$country2==i,], model="random")
  s <- summary(model)$coef
  coefs$b[coefs$countries==i] <- s[9,1]
  coefs$se[coefs$countries==i] <- s[9,2]
}
coefs$low <- coefs$b - (coefs$se * 1.96)
coefs$high <- coefs$b + (coefs$se * 1.96)
coefs$range <- as.numeric(rownames(coefs))

{
pdf("figure3.pdf", width = 10, height = 5, bg="white") 
par(mar=c(5.1,4.1,1.1,2.1))
plot(coefs$range,coefs$b,type="n",ann=F,axes=F,ylim=c(-1.5,1.5))
abline(h=0,lwd=1,lty=2,col="black")
arrows(coefs$range,coefs$low,coefs$range,coefs$high,angle=90,length=.05,code=3,col="gray",lwd=2,lty=1)
points(coefs$range,coefs$b,pch=15,cex=1.5)
axis(1,col="white",at=c(1:19),labels=coefs$countries,cex.axis=.8)
axis(2,las=1)
mtext("Ideological Difference in the Effect of Negativity",side=2,line=3)
invisible(dev.off())
}

```

```{r STATA stimulus-data analyses, as.is=T}

stata_src <- readLines("Stimulus-data.do")
stata(stata_src,data.in=R)

```

## Analysis of Time-series-data-videos.dta

```{r Figure 4}

coefs <- as.data.frame(countries)
coefs$b <- NA
coefs$se <- NA

#V$neg <- V$negativity
A <- V[V$period==2 & V$local==0,]

Pd <- pdata.frame(A,index = c("resp", "timesec"))
Pd$Lgslmc <- lag(Pd$gslmc,1)
Pd$Dgslmc <- Pd$gslmc - Pd$Lgslmc
for (i in countries){
model <- plm(Dgslmc ~ negativity * timesecStorylog + negativity * wp_right2  + Lgslmc 
             + female + age + income + university + order,
              data=Pd[Pd$country2==i,], model="random") 
  s <- summary(model)$coef
  coefs$b[coefs$countries==i] <- s[12,1]
  coefs$se[coefs$countries==i] <- s[12,2]
}

coefs$low <- coefs$b - (coefs$se * 1.96)
coefs$high <- coefs$b + (coefs$se * 1.96)
coefs$range <- as.numeric(rownames(coefs))

{
pdf("figure4.pdf", width = 10, height = 5, bg="white") 
par(mar=c(5.1,4.1,1.1,2.1))
plot(coefs$range,coefs$b,type="n",ann=F,axes=F,ylim=c(-.025,.025))
abline(h=0,lwd=1,lty=2,col="black")
arrows(coefs$range,coefs$low,coefs$range,coefs$high,angle=90,length=.05,code=3,col="gray",lwd=2,lty=1)
points(coefs$range,coefs$b,pch=15,cex=1.5)
axis(1,col="white",at=c(1:19),labels=coefs$countries,cex.axis=.8)
axis(2,las=1)
mtext("Ideological Difference in the Effect of Negativity",side=2,line=3.2)
invisible(dev.off())
}

```

```{r STATA time-series-data-videos analyses, as.is=T}

stata_src <- readLines("Time-series-data-videos.do")
stata(stata_src,data.in=V)

```

## Analysis of Time-series-data-photos.dta

```{r STATA time-series-data-photos analyses, as.is=T}

P$country2 <- as.character(P$country2)
stata_src <- readLines("Time-series-data-photos.do")
stata(stata_src,data.in=P)

```



